home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_3.7 / multires / multires.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  14.0 KB  |  383 lines
  1. /* 
  2.  * multires.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* MULTIRES:    program performs image subsampling after low-pass filtering
  10.  *            to produce multiresolution pyramid of images subsampled
  11.  *              by 2. These are displayed in the original image area.
  12.  *                  usage: multires inimg outimg [-c || -m || -l LEVEL]
  13.  *                                    [-d DEGREE_FILTERING] [-q quickFlag] [-L]
  14.  */
  15.  
  16. #define OUTTYPE_DFLT 1          /* =1 composite image; =2 multiple; =3 single level */
  17. #define BGFLAG_DFLT 0           /* background flag for composite =0 blank; =1 orig. */
  18. #define SUBRATE 2               /* subsample rate between levels */
  19. #define DFLT_DEG_FLTR 0.5       /* default degree of filtering */
  20. #define RECT_CUTOFF 0.443883    /* factor to calc. 3dB cutoff freq */
  21. #define GAUSS_CUTOFF 0.832555   /* factor to calc. 3dB cutoff freq */
  22. #define ROUNDUP 0.99999         /* for rounding up integer */
  23. #define NORM 1000.0             /* normaliz. factor to integerize */
  24.  
  25. #include <stdio.h>
  26. #include <stdlib.h>
  27. #include <images.h>
  28. #include <tiffimage.h>
  29. #include <math.h>
  30. #include <string.h>
  31. extern void print_sos_lic ();
  32.  
  33. long convolve1r (Image *, Image *, long *, long, long);
  34. long smoothr (Image *, Image *, long, long);
  35. long input (int, char **, short *, long *, short *, double *, short *);
  36. long usage (short);
  37.  
  38. main (argc, argv)
  39.      int argc;
  40.      char *argv[];
  41. {
  42.   Image *imgI, *imgO;           /* input, output image structures */
  43.   Image *imgT;                  /* intermediate image structure */
  44.   unsigned char **imgIn, **imgOut;  /* input, output image arrays */
  45.   unsigned char **imgTi;        /* intermediate image array */
  46.   char outFile[80];             /* output filenames, multipe levels */
  47.   long width, height;           /* image size */
  48.   double cutoff;                /* low-pass cutoff frequency */
  49.   long fltrSize;                /* filter mask size */
  50.   short outType;                /* =1 composite image; =2 multiple; =3 single level */
  51.   long outLevel;                /* for single level output, specifies which level */
  52.   short bgFlag;                 /* background flag for composite =0 blank; =1 orig. */
  53.   double degFltr;               /* degree of filtering [0.0 - 1.0] */
  54.   double stdDev;                /* std. deviation param. of Gaussian */
  55.   long *fltr;                   /* 1D filter mask */
  56.   long midFltr;                 /* middle filter coefficient */
  57.   long widthR, heightR;         /* subsampled image sizes */
  58.   short quickFlag;              /* for quick and dirty rect. filter */
  59.   long xOffSet;                 /* x offset for pyramid display */
  60.   long x, xO, y, i;
  61.  
  62.   if (input (argc, argv, &outType, &outLevel, &bgFlag,
  63.              °Fltr, &quickFlag) < 0)
  64.     return (-1);
  65.  
  66. /* open input image */
  67.   imgI = ImageIn (argv[1]);
  68.   imgIn = ImageGetPtr (imgI);
  69.   height = ImageGetHeight (imgI);
  70.   width = ImageGetWidth (imgI);
  71.   printf ("Image size is %dx%d.\n", width, height);
  72.  
  73. /* determine filter size and coefficients */
  74.   if (degFltr == 0.0)
  75.     cutoff = 100.0;
  76.   else
  77.     cutoff = 0.5 / (SUBRATE * degFltr);
  78.   if (quickFlag)                /* rectangular filter window */
  79.     fltrSize = (long) (RECT_CUTOFF / cutoff + ROUNDUP);
  80.   else {                        /* Gaussian filter window */
  81.     stdDev = GAUSS_CUTOFF / cutoff;
  82.     fltrSize = (long) (3.0 / (2.0 * cutoff) + 0.5);
  83.   }
  84.   if ((fltrSize % 2) == 0)
  85.     fltrSize += 1;
  86.   printf ("Filter 3dB cutoff frequency is %5.2f and filter size is %dx%d.\n",
  87.           cutoff, fltrSize, fltrSize);
  88.  
  89.   if (!quickFlag) {
  90.     fltr = (long *) calloc (fltrSize, sizeof (*fltr));
  91.     midFltr = fltrSize / 2;
  92.     for (i = 0; i <= midFltr; i++)
  93.       fltr[i + midFltr] = fltr[midFltr - i]
  94.         = (long) (NORM * exp (-i * i / (2.0 * stdDev * stdDev)));
  95.   }
  96.  
  97. /* construct multi-resolution pyramid */
  98.   heightR = height;
  99.   widthR = width;
  100.   imgT = ImageAlloc (height, width, 8);
  101.   imgTi = ImageGetPtr (imgT);
  102.   imgO = ImageAlloc (heightR, widthR, 8);
  103.   imgOut = ImageGetPtr (imgO);
  104.   xOffSet = 0;
  105.   for (y = 0; y < height; y++)
  106.     for (x = 0; x < width; x++)
  107.       imgTi[y][x] = imgIn[y][x];
  108.   if (outType == 1 && bgFlag == 0) {
  109.     for (y = 0; y < height; y++)
  110.       for (x = 0; x < width; x++)
  111.         imgIn[y][x] = 0;
  112.   }
  113.   for (i = 1; (heightR >= 4 && widthR >= 4); i++) {
  114.     if (quickFlag)
  115.       smoothr (imgT, imgO, fltrSize, SUBRATE);
  116.     else
  117.       convolve1r (imgT, imgO, fltr, fltrSize, SUBRATE);
  118.     heightR = heightR / SUBRATE;
  119.     widthR = widthR / SUBRATE;
  120.     printf ("level %d: %dx%d\n", i, heightR, widthR);
  121.     ImageFree (imgT);
  122.     imgT = ImageAlloc (heightR, widthR, 8);
  123.     imgTi = ImageGetPtr (imgT);
  124.     for (y = 0; y < heightR; y++)
  125.       for (x = 0, xO = xOffSet; x < widthR; x++, xO++)
  126.         imgIn[y][xO] = imgTi[y][x] = imgOut[y][x];
  127.     if (outType == 2) {         /* multiple level output images */
  128.       sprintf (outFile, "m%d_%s", i, argv[2]);
  129.       printf ("filename is %s\n", outFile);
  130.       ImageOut (outFile, imgT); /* this is a hack -- flipBits in */
  131. //  ImageOut (outFile, imgT); /* tiffimage.c causes inversion */
  132.     }
  133.     else if (i == outLevel) {   /* single level output image */
  134.       ImageOut (argv[2], imgT);
  135.       exit (1);
  136.     }
  137.     xOffSet += widthR;
  138.   }
  139.  
  140. /* output subsampled image */
  141.   ImageOut (argv[2], imgI);
  142.  
  143.   return (0);
  144. }
  145.  
  146.  
  147.  
  148. /* SMOOTHR:     function performs 2-D smoothing by separable, 1-D averaging
  149.  *            and returns subsampled image
  150.  *                  usage: smoothr (imgI, imgO, nFltr, rSubsample)
  151.  */
  152.  
  153. long
  154. smoothr (imgI, imgO, nFltr, rSubsample)
  155.      Image *imgI, *imgO;        /* input,output image structures */
  156.      long nFltr;                /* number of filter coefficients */
  157.      long rSubsample;           /* subsample rate */
  158. {
  159.   Image *imgT;                  /* intermediate image structure */
  160.   unsigned char **imgIn,        /* input image */
  161.   **imgTi,                      /* intermediate image */
  162.   **imgOut;                     /* output image */
  163.   long width, height;           /* size of image */
  164.   long sum;                     /* sum of filter convolution at a pixel */
  165.   long midFltr;                 /* middle coefficient index of filter */
  166.   long xEnd, yEnd;              /* end coefficients of convolution */
  167.   long xOut, yOut;              /* output image coordinates */
  168.   long x, y, i;
  169.  
  170. /* initialize images */
  171.   imgIn = ImageGetPtr (imgI);   /* input image */
  172.   height = ImageGetHeight (imgI);
  173.   width = ImageGetWidth (imgI);
  174.   imgT = ImageAlloc (height, width, 8);  /* intermediate image */
  175.   imgTi = ImageGetPtr (imgT);
  176.   imgOut = ImageGetPtr (imgO);  /* output image */
  177.  
  178. /* perform row-wise smoothing */
  179.   midFltr = (nFltr - 1) / 2;
  180.   xEnd = width - midFltr;
  181.   for (y = 0; y < height; y++) {
  182.     for (x = midFltr; x < xEnd; x += rSubsample) {
  183.       sum = imgIn[y][x];
  184.       for (i = 1; i <= midFltr; i++)
  185.         sum += (long) (imgIn[y][x - i] + imgIn[y][x + i]);
  186.       imgTi[y][x] = (unsigned char) (sum / nFltr);
  187.     }
  188.   }
  189.  
  190. /* perform column-wise convolution */
  191.   yEnd = height - midFltr;
  192.   for (y = midFltr, yOut = 0; y < yEnd; y += rSubsample) {
  193.     for (x = midFltr, xOut = 0; x < xEnd; x += rSubsample) {
  194.       sum = imgTi[y][x];
  195.       for (i = 1; i <= midFltr; i++)
  196.         sum += imgTi[y - i][x] + imgTi[y + i][x];
  197.       imgOut[yOut][xOut++] = (unsigned char) (sum / nFltr);
  198.     }
  199.     yOut++;
  200.   }
  201.  
  202.   return (0);
  203. }
  204.  
  205.  
  206.  
  207. /* CONVOLVE1R:  function performs 2-D convolution by separable, 
  208.  *            symmetric, 1-D filter mask, then outputs subsample image 
  209.  *                  usage: convolve1r (imgI, imgO, fltr1D, nFltr)
  210.  */
  211.  
  212. long
  213. convolve1r (imgI, imgO, fltr1D, nFltr, rSubsample)
  214.      Image *imgI;               /* input image structure */
  215.      Image *imgO;               /* output image structure */
  216.      long *fltr1D;              /* 1-D of filter mask */
  217.      long nFltr;                /* number of filter coefficients */
  218.      long rSubsample;           /* subsample rate */
  219. {
  220.   Image *imgT;                  /* intermediate image structure */
  221.   unsigned char **imgIn,        /* input image */
  222.   **imgTi,                      /* intermediate image */
  223.   **imgOut;                     /* output image */
  224.   long width, height;           /* size of image */
  225.   long sumFltr;                 /* sum of fltr coefficients */
  226.   long sum;                     /* sum of filter convolution at a pixel */
  227.   long midFltr;                 /* middle coefficient index of filter */
  228.   long xEnd, yEnd;              /* end coefficients of convolution */
  229.   long xOut, yOut;              /* output image coordinates */
  230.   long x, y, i;
  231.  
  232. /* initialize images */
  233.   imgIn = ImageGetPtr (imgI);
  234.   height = ImageGetHeight (imgI);
  235.   width = ImageGetWidth (imgI);
  236.   imgT = ImageAlloc (height, width, 8);
  237.   imgTi = ImageGetPtr (imgT);
  238.   imgOut = ImageGetPtr (imgO);
  239.  
  240. /* find sum of filter coefficients */
  241.   for (i = 0, sumFltr = 0; i < nFltr; i++)
  242.     sumFltr += fltr1D[i];
  243.  
  244. /* perform row-wise convolution */
  245.   midFltr = (nFltr - 1) / 2;
  246.   xEnd = width - midFltr;
  247.   for (y = 0; y < height; y++) {
  248.     for (x = midFltr; x < xEnd; x += rSubsample) {
  249.       sum = fltr1D[midFltr] * imgIn[y][x];
  250.       for (i = 1; i <= midFltr; i++)
  251.         sum += fltr1D[i + midFltr]
  252.           * (imgIn[y][x - i] + imgIn[y][x + i]);
  253.       imgTi[y][x] = (unsigned char) (sum / sumFltr);
  254.     }
  255.   }
  256.  
  257. /* perform column-wise convolution */
  258.   yEnd = height - midFltr;
  259.   for (y = midFltr, yOut = 0; y < yEnd; y += rSubsample) {
  260.     for (x = midFltr, xOut = 0; x < xEnd; x += rSubsample) {
  261.       sum = fltr1D[midFltr] * imgTi[y][x];
  262.       for (i = 1; i <= midFltr; i++)
  263.         sum += fltr1D[i + midFltr] * (imgTi[y - i][x] + imgTi[y + i][x]);
  264.       imgOut[yOut][xOut++] = (unsigned char) (sum / sumFltr);
  265.     }
  266.     yOut++;
  267.   }
  268.  
  269.   return (0);
  270. }
  271.  
  272.  
  273.  
  274. /* USAGE:       function gives instructions on usage of program
  275.  *                    usage: usage (flag)
  276.  *              When flag is 1, the long message is given, 0 gives short.
  277.  */
  278.  
  279. long
  280. usage (flag)
  281.      short flag;                /* flag =1 for long message; =0 for short message */
  282. {
  283.  
  284.   printf ("USAGE: multires inimg outimg [<-c> || <-m> || <-l LEVEL>]\n");
  285.   printf ("                             [-b] [-d DEG_FILTERING <%3.1f>] [-q] [-L]\n", DFLT_DEG_FLTR);
  286.   if (flag == 0)
  287.     return (-1);
  288.  
  289.   printf ("\nmultires performs multi-resolution processing upon an image\n");
  290.   printf ("ARGUMENTS:\n");
  291.   printf ("    inimg: input image filename (TIF)\n");
  292.   printf ("   outimg: output image filename (TIF)\n\n");
  293.   printf ("OPTIONS:\n");
  294.   printf ("                -c: composite image output of all multi-resolution levels\n");
  295.   printf ("                    arranged in a single image; this is the default mode.\n");
  296.   printf ("                -m: multiple images of each multi-resolution level\n");
  297.   printf ("                    in separate output image files; the output files are\n");
  298.   printf ("                    m1_outimg.tif, m2_outimg.tif, ... for each level 1, 2, etc.\n");
  299.   printf ("          -l LEVEL: output is single image of the multi-resolution\n");
  300.   printf ("                    level specified; level 1 is half size of original, level 2\n");
  301.   printf ("                    is half size of level 1, etc; filename is output\n");
  302.   printf ("                    filename specified, with 1, 2, etc. appended.\n");
  303.   printf ("                    NOTE: only one of -c, -m, or -l should be chosen.\n");
  304.   printf ("                -b: BACKGROUND_FLAG for composite image <-c>; if flag is\n");
  305.   printf ("                    set, then background to composite images is original image;\n");
  306.   printf ("                    otherwise, the default is a blank background\n");
  307.   printf ("                    of the original image size.\n");
  308.   printf ("  -d DEG_FILTERING: value set between 0.0 and 1.0 controlling\n");
  309.   printf ("                    the degree of low-pass filtering; a value\n");
  310.   printf ("                    of 1.0 produces maximum filtering to reduce aliasing,\n");
  311.   printf ("                    but may lead to blurring;\n");
  312.   printf ("                    a value of 0.0 does no filtering at the risk of aliasing.\n");
  313.   printf ("                    the default is %3.1f.\n", DFLT_DEG_FLTR);
  314.   printf ("                -q: when set, performs quick low-pass filtering \n");
  315.   printf ("                    with rectangular filter rather than Gaussian filter;\n");
  316.   printf ("                    the former is faster but doesn't produce equally\n");
  317.   printf ("                    good results as does the latter.\n");
  318.   printf ("                -L: print Software License for this module\n");
  319.  
  320.  
  321.   return (-1);
  322. }
  323.  
  324.  
  325. /* INPUT:       function reads input parameters
  326.  *               usage: input (argc, argv, &outType, &outLevel, &bgFlag,
  327.  *                                              °Fltr, &quickFlag)
  328.  */
  329.  
  330. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  331.  
  332. long
  333. input (argc, argv, outType, outLevel, bgFlag, degFltr, quickFlag)
  334.      int argc;
  335.      char *argv[];
  336.      short *outType;            /* =1 composite image; =2 multiple; =3 single level */
  337.      long *outLevel;            /* for single level output, specifies which level */
  338.      short *bgFlag;             /* background flag for composite =0 blank; =1 orig. */
  339.      double *degFltr;           /* degree of filtering [0.0 - 1.0] */
  340.      short *quickFlag;          /* for quick and dirty rect. filter */
  341. {
  342.   long n;
  343.  
  344.   if (argc < 3)
  345.     USAGE_EXIT (1);
  346.  
  347.   *outType = OUTTYPE_DFLT;
  348.   *outLevel = 0;
  349.   *bgFlag = BGFLAG_DFLT;
  350.   *degFltr = DFLT_DEG_FLTR;
  351.   *quickFlag = 0;
  352.  
  353.   for (n = 3; n < argc; n++) {
  354.     if (strcmp (argv[n], "-c") == 0)
  355.       *outType = 1;
  356.     else if (strcmp (argv[n], "-m") == 0)
  357.       *outType = 2;
  358.     else if (strcmp (argv[n], "-l") == 0) {
  359.       if (++n == argc || argv[n][0] == '-')
  360.         USAGE_EXIT (0);
  361.       *outLevel = atol (argv[n]);
  362.       *outType = 3;
  363.     }
  364.     else if (strcmp (argv[n], "-b") == 0)
  365.       *bgFlag = 1;
  366.     else if (strcmp (argv[n], "-d") == 0) {
  367.       if (++n == argc || argv[n][0] == '-')
  368.         USAGE_EXIT (0);
  369.       *degFltr = atof (argv[n]);
  370.     }
  371.     else if (strcmp (argv[n], "-q") == 0)
  372.       *quickFlag = 1;
  373.     else if (strcmp (argv[n], "-L") == 0) {
  374.       print_sos_lic ();
  375.       exit (0);
  376.     }
  377.     else
  378.       USAGE_EXIT (0);
  379.   }
  380.  
  381.   return (0);
  382. }
  383.